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ABSTRACT 

A considerable fraction of multi-planet systems discovered by the observational surveys of extra- 
solar planets reside in mild proximity to first-order mean motion resonances. However, the relative 
remoteness of such systems from nominal resonant period ratios (e.g. 2:1, 3:2, 4:3) has been inter- 
preted as evidence for lack of resonant interactions. Here we show that a slow divergence away from 
exact commensurability is a natural outcome of dissipative evolution and demonstrate that libration 
of critical angles can be maintained tens of percent away from nominal resonance. We construct an 
analytical theory for the long-term dynamical evolution of dissipated resonant planetary pairs and 
confirm our calculations numerically. Collectively, our results suggest that a significant fraction of the 
near-commensurate extrasolar planets are in fact resonant and have undergone significant dissipative 
evolution. 



1. INTRODUCTION 

Among the most unexpected discoveries brought forth 
by extrasolar planetary surveys to date has been the 
identification of numerous planetary bodies that reside 
in close proximity to their host stars. Planets of this sort 
are of great scientific interest because they represent a 
class of objects unavailable for study in our own solar 
system. In turn, observational characterization of such 
planetary systems can yield avenues towards identifying 
specific physical/dynamical behavior that does not occur 
locally, thus broadening our knowledge of the possible 
evolutions of planetary systems. 

A readily apparent dynamical feature of close-in extra- 
solar planetary systems, highlighte d by observational 
surveys such as the K epler mission ( |Borucki et alpoTT) 



Lissauer et al. 2011), is the prominence of near mean- 
motion commensurabilities (i.e. integer period ratios) 
among sub-giant planets (Figure 1). Accordingly, un- 
derstanding how close-in planetary systems attain near- 
resonant orbital architectures is the primary focus of this 
work. 

The process of resonant locking requires slow, con- 



vergent orbital ev olution of planetary bodies (Goldreich 



1965| |Peale 19761). It is likely that torques associated 
with disk-driven migration often lead to resonant cou- 
pling, and it has been suggested that near-exact com- 
mensurability should be maintained as planets travel 
through their proto-planetary disks ([Terquem fc Pa- 
paloizou||2007l iCresswell & Nelson||2008J) . However, the 



onset of magneto-rotational instability ( fBalbus fc Haw- 
|ley||199ij ) and the associated turbulence in protoplane- 
tary disks can act to disrupt mean-motion resonances 
d Adams et al.||2008[ IRein & Papaloizou||2009| IKetchum 



et al.||2011 



hus, if disks are violently turbulent, reso 
nant objects should be rare. 

As already hinted above, the observations show that 
there exists a characteristic regime in between the two ex- 
tremes, and the precise dynamical nature of this regime 
is elusive. Particularly, planets often reside sufficiently 
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far away (a few percent or more) from their nominal 
first-order resonant locations (i.e. period ratios of 2:1, 
3:2, 4:3) to be readily interpreted as non-resonant. Yet 
the preference for orbits just wide of resonance and a 
characteristic pile- up of near-resonant objects (Fig. 1) is 
suggestive of a common evolutionary path. Indeed, the 
mechanism responsible for such configurations has been 



noted to be a sub ject of great theoretical interest (Fab 



rycky et al. 2012) 



It is possible in principle that most sub-giant planets 
arrive onto their close-in orbits in resonance and subse- 
quently diverge away from exact commensurability due 
to tidal dissipation. Tides alone affect the semi-major 
axes only on very long timescales (often much longer 
than the Hubble time). However, as shown by the non- 
linear perturbative calculations and A-body simulations 
aimed at reproducing the orbital con figurations of the 
HD40307 dPapaloizou fc Terquem|20l0l ) as well as GL581 
and HD10180 (Papaloizou 2011 ) systems, resonant in- 
teractions can be quite effective at converting tidal ec- 
centricity damping (which acts much faster) into a di- 
vergence of the orbital semi-major axes of the resonant 
bodies. In particular, the said simulations suggest that 
resonant coupling can be maintained far from nominal 
resonant locations and significantly aids in enhancing or- 
bital divergence. 
The calculations performed by ( |Papaloizou fc Terquem| 



2010 1 motivate our development of a general qualitative 
understanding of the orbital evolution of close-in res- 
onant planetary systems subject to dissipative effects. 
Thus, the development of an analytical theory for dissi- 
pative divergence of resonant orbits is the primary focus 
of this paper. The number of well-characterized systems 
within the Kepler sample remains limited and estimation 
of planetary masses from radii alon e is generally risky 
( Stevensoii]|1982[ |Rogers et al. 2011). Consequently, in 
this work, we shall concentrate our efforts on characteri- 
zation of the physical process rather than reproduction of 
any particular orbital architecture. Still, we argue that 
the interplay between resonant effects and tidal dissipa- 
tion is the primary mechanism by which planets attain 
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nominal resonance location 
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Fig. 1. — A histogram of the period ratios of all planet pairs de- 
tected by the Kepler mission with no niters on planetary radius 
or orbital period (http://planetquest.jpl.nasa.gov/kepler). In sys- 
tems where more than two planets are present, only the neighboring 
period ratios are reported. Note the highlighted enhancement of 
objects immediately outside of the common (2:1 and 3:2) first-order 
mean motion resonances. 



near-commensurate orbits. Lithwick & Wu (2012) ar- 
rived at many of the results presented m this work si- 
multaneously and independently; their paper was posted 
on arxiv.org at the same time as this one. 

The paper is organized as follows. In section 2, we set 
the stage by developing an integrable approximation to 
the conservative dynamics of a resonant pair at low ec- 
centricities and validate the theory by comparison with 
A^-body simulations. In section 3, we introduce dissipa- 
tion into the problem and show that tidal effects drive 
the system towards a quasi-stationary state that is char- 
acterized by an irreversible drift away from nominal res- 
onance, where the inner planet's orbit decays at a rate 
that is faster than that expected from the direct tidal 
effect, while the outer planet gains orbital energy. In 
section 4, we discuss the extension of our formalism to 
multi-resonant systems. Subsequently, we conclude and 
discuss our results in section 5. 

2. CONSERVATIVE DYNAMICS OF A RESONANT 
PLANETARY PAIR 

Resonant dynamics of planetary pairs have be en stud- 
ied by numerous a uthors in the past (see Ch.8 of|Murray 



& Dermott ( 1999 ) and the references therein) . This work 
builds on their contributions. 

Our eventual goal is to construct an analytical model 
for the long-term evolution of resonant orbits under dis- 
sipative effects. Before complicating the picture with dis- 
sipation, however, we must first build a purely analytical 
model for conservative resonant interactions. Thus, in 
this section, we shall derive a simple, physically intu- 
itive closed-form solution for the time-evolution of a res- 
onant planetary pair. Accordingly, we shall firs t work in 
the sp irit of classical perturbation theory (e.g. [Message] 
( 1966); [Peale| (1986)) and employ numerical calculations 
primarily as a means of confirmation. 

Let us begin by considering a quasi-intcgrable Hamil- 
tonian of the form 



where 



H = H iep + H res + 0{e 2 ,n 
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is the Keplerian Hamiltonian and 
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is the first-order k : k — 1, k G Z resonant perturbation. 
Here, the orbital elements take on their standard nota- 
tion, M is the mass of the central star and TOi,m 2 are 
the masses of the planets with the subscript 1 and 2 re- 
ferring to the inner and outer planets respectively. The 

quantities /res and / r e S depend on the semi-major axis 
ratio (ai/a 2 ) only and are tabulated in the literature (see 
for example Murray & Dermott (1999)). 

Because Keplerian orbital elements are not canonically 
conjugated, we revert to Poincare variables for further 
calculations: 



A = mVGMa, 
T = A(1- y/l- 
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e 2 ) « A e 2 /2, 7 



-VJ, 
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where Af is the mean anomaly and the indexe 1, 2 are 
omitted for simplicity. In terms of the Poincare variables, 
the Hamiltonians, "Hkcp and T-L res read: 
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As already implied by equation ([!]), we shall work to 
first order in eccentricity, neglecting secular effects and 
resonances of order greater than unity. Generally, H only 
constitutes a good approximation to the true dynamics 
of a planetary pair in the vicinity of a mean-motion res- 
onance. 

Because the perturbation "H ros is of order e, we expect 
that the semi major axes can change by 0(y/e) relative 
to their nominal, resonant values. Thus, we expand the 
terms in "Hkep to second order in SA = A — [A] , where [A] 
is the nominal value of A: 
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Consistently, we evaluate T~L rcs in (6) at [A], as it is al- 
ready of order 0(e). Constant terms are dynamically 
unimportant and can thus be dropped from the Hamil- 
tonian, implying 5 A A and <5A 2 — !> A 2 - 2 A [A]: 
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Fig. 2. — Orbital evolution of a nearly mass-less (m = 10 10 Mq) particle in an interior 2:1 mean motion resonance with a Jupieter-mass 
object (m = 10 — 3 Mq) with a semi-major axis of a = 1 AU. The evolution is shown over 50 orbital periods of the perturbing object, 
corresponding to ~ 5 resonant cycles. The panels A, B, and C show the variation in the particle's semi-major axes, eccentricity and the 
critical resonant angle respectively. The red curve was obtained analytically utilizing the framework developed in section 2. The blue curve 
was obtained by numerically integrating the equations of motion that arise from the Hamiltonians j(3]l and Q. The gray curve is a result 
of a direct N-body simulation. 



Note that the planetary mean motion is given by 



dX 



dU 
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As a result, "Hkep can be rewritten in a compact form: 

H kcp = 4([n]iAi + [n] 2 A a ) - |([fc]iA? + [ft] 2 A|), (11) 

where [h] = [n]/[A] = m/[a] 2 . 

Although 'Hkcp is now expressed in a simple form, H ies 
remains cumbersome largely due to the formulation of 
the resonant angles which appear as cosine arguments. 
Let us employ a canonical transformation of coordinates, 
utilizing the following generating function of the second 
kind: 

F 2 = Ai*i + A 2 * 2 + (fcA 2 — (k - l)Ai + 7i)$i 
+ (fcA 2 -(fc-l)A 1 + 72 )$ 2 , (12) 

where and $ are new momenta. Upon application of 
the transformation equations 

OF dF 
dA ~~ &i 



A 



(13) 



we obtain new canonically conjugated action-angle vari- 
ables 

*i =Ai + (k - l)($ a + $ 2 ) ^i=Ai 
* 2 =A 2 -A;($i + $2) 1P2 = A 2 

$i =Ti <j) X = k\ 2 - (k - 1)A X + 71 
$ 2 =r 2 ^ 2 = fcA 2 - (A - l)Ai + 72. (14) 

In terms of these variables, the resonant contribution to 
H is expressed as follows: 
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while the Keplerian contribution reads: 
H kcp =4[n] 1 (* 1 -(fc-l)($ 1 + 
+ 4[n] 2 (* 2 + fc($i +$2)) 



* 2 )) 
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The transformation to new variables allows us to make 
further simplifications to %kep- Specifically, because 
dH/dtp = 0, ^1 and iff 2 are constants of motion, allow- 
ing us to drop additional terms. It is further instructive 
to recall that 4> oc e 2 . Consequently, if e -C 1, non- 
linear terms proportional to $ 2 , $21 an d ^1^2 can be ne- 
glected. This approximation filters out chaotic dynamics 
from the Hamiltonian and therefore will not yield an ad- 
equate representation of the ev olution of the system in 



the r esonances overlap region (Chirikov 1979 
19801. However as will be shown below. 



Wisdom 
this assump 

tion is well satisfied in the calculations of interest. Upon 
making these simplifications, the Keplerian Hamiltonian 
is simply 

Wte P = (4(fc[n] 2 -(fc-l)[n] 1 )) 

+ 3{[h]i(k - 1)*! - [A] 2 A* 2 ))(*i + $2). (17) 

Note that by definition, (k[n]2 — (k— l)[n]i) = because 
it signifies exact resonance. As a result, only terms pro- 
portional to [h] remain in "Hkcp- 
The full Hamiltonian now takes on a very simple form: 



where 



r) = 3([h] 1 (k-l)y 1 -[h] 2 M 2 ) 



(18) 
(19) 

is related the circulation frequency of the critical angles 
in an unperturbed case (mi = m,2 — 0) and is thus a 
measure of proximity of the planetary pair to exact Kep- 
lerian resonance (note that rj — > as A — > [A] and $ — » 0, 
corresponding to if> = [A]) while 
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(20) 



are the strengths of t he resonances. It is noteworthy that 
the Hamiltonian (18) represents two decoupled Hamilto- 
nians, each of which has a form similar of the "second 
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fundamental model of resonance" (Henrard & Lamaitre 



19831, apart from the missing term, proportional to 4> z 
that we have neglected. 

In the coordinates used up to now, the equations of 
motion are singular at $ = 0. However, this singularity 
can be overcome by switching to mixed cartesian coordi- 
nates 

x = %/2$sin(» y = V2$cos(<f>) (21) 

via a contact transformation (here, x is identified as the 
coordinate and y as the momentum). In these coordi- 
nates, the Hamiltonian reads 
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Accordingly, the equations of motion are: 
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Fig. 3. — Orbital evolution of a nearly mass-less (m = 10 10 Mq) particle in an exterior 3:2 mean motion resonance with a Jupieter-mass 
object (m = 10 — 3 Mq) with a semi-major axis of a = 1 AU. The evolution is shown over 50 orbital periods of the perturbing object, 
corresponding to ~ 7 resonant cycles. As in Figure (2), the panels A, B, and C show the variation in the particle's semi-major axes, 
eccentricity and the critical resonant angle respectively. The red curve was obtained analytically utilizing the framework developed in 
section 2. The blue curve was obtained by numerically integrating the equations of motion that arise from the Hamiltonians (JgJ and 
The gray curve is a result of a direct N-body simulation. 

where C\ and C2 are (possibly complex) constants of inte- 
gration. Note that except for a dependence of the leading 
term on z, equations (25) are analogous to the complex 
formulation of the Laplace-Lagrange theory for secular 
intera ctions (Wu & Goldrcich 20 02| |Batygin fc Laughlin 
2011 1, although the variables take on a different meaning. 

Within the context of this model, variations in semi- 
major axes can be derived from the fact that ^> remain 
constants of motion. Examples of the application of the 
theory are presented in Figures (2) and (3). In both of 
the illustrated cases, a nearly mass-less (m = 1O~ 1O M0) 
particle is perturbed by a Jupieter-mass object (m — 
1O _3 M0) with a semi-major axis of a = 1 AU. Figure 
(2) shows an interior 2:1 mean motion resonance while 
Figure (3) shows an exterior 3:2 mean motion resonance. 
The red curves denote analytical theory, the blue curves 
represent a numerical integration of the non-linear per- 
turbative Hamiltonians ^ and (FF|), and the gray curves 
are the results of numerical A-Dody simulations, per- 
formed using the hybrid algorith m of the orbital in tegra- 
tion software package mercury6 ( |Chambers|[l99 9 ) . Note 
that as a consequence of the simplifications made in or- 
der to express the analytical solution in closed form, the 
blue (non-linear perturbative) curve has slightly differ- 
ent frequency and amplitude of oscillation relative to the 
red (analytical) curve, although the two curves exhibit 
the same qualitative behavior. However, in addition to 
the resonant variations, the grey (A-body) curve shows 
non-resonant, short-period oscillations, that are filtered 
out by retaining only the resonant terms in the Hamilto- 
nian. These short-periodic oscillations are unimportant 
to the problem at hand, as they do not contribute to the 
time- averages of the resonant angles. Note also that, al- 
though the particles in both examples are relatively far 
away from nominal resonance, the critical angles remain 
in libration. 

3. DISSIPATIVE DYNAMICS OF A RESONANT 
PLANETARY PAIR 

There exists an abundance of circumstances where the 
evolution of a planetary system cannot be described in 
terms of strictly conservative interactions. For exam- 
ple, planets embedded in protoplanetary disks exp erience 
dissipative forces exerted by the gaseous nebula (Lee & 
Peale|[2002 ), while planets that reside on orbits that are 
in close proximity to their host stars are subject to tidal 
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Although we can continue to work in terms of the mixed 
cartesian coordinates, the equations of motion can be 
re-written in a more compact form by treating x and y 
as imaginary and real components of a single complex 
variable 

z — ix + y. (24) 

Now, the equations of motion can be written down con- 
cisely: 

dzi 

' H + ITjZi 
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Fig. 4. — Dissipative evolution of an equal-mass (mi = = 10~ 4 M) planetary pair in a 2:1 mean motion resonance over t/r = 100 
circularization timescales. Panels A and D show the evolution of the planetary semi-major axes. Note that at all times dissipative 
interactions give rise to a monotonic divergence of the orbits. This can be further inferred from panel E which shows the measure of 
proximity to exact resonance, r) < monotonically decreasing. Panels C and F show the evolution of the critical angles. Note that the 
system attains a state of quasi-cquilibrium after f ~ 5t. Accordingly, the eccentricity evolution also becomes quasi-stationary after the 
critical angles collapse to a near-focal state. The red curves were obtained analytically utilizing the framework developed in section 3. The 
blue curves were obtained by numerically integrating the equati ons of mo tion that arise from the Hamiltonians and J?), augmented 
with a simple parameterization of tidal dissipation (i.e. equations \27\ and |3T}). The gray curv es were computed num erically v 
N-body simulation where dissipation has been taken into account using the tidal framework of | |Eggleton et al.|1998[ l. 



ify with a direct 



friction ( |Bodenheimer et al.|2001 ) (in this work, we shall 
concentrate on the latter). In the extrasolar context, 
tidal dissipation usually results in the decay of orbital 
eccentricity and semi-major axes. 

With the exception of special configurations, the char- 
acteristic timescales for the decay of eccentricity and 
semi-major axes differ significantly (often by orders of 
magnitude). This is in part because the changes in ec- 
centricity are controlled by the rate of angular momen- 
tum exchange in the system, while changes in the semi- 
major axes are largely governed by the rate of energy 
dissipation, which is usually a much slower process. As 
a result for the purposes of this work, we shall invoke 
separation of timescales and treat the decays of e and 
a independently. For e <C 1, the or bit-averaged rate of 
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tidal eccentricity decay is given by ( |Goldreich fe Soter 
19661): 
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where k is the planetary Love number, Q is the tidal 
quality factor (note that dissipation within the host-star 
is neglected as usual), and R is the planetary radius. 
Noting that \z\ ~ e-v/[A], it is trivial to incorporate ec- 
centricity decay into equations ( 25 ) : 

dzi 

' VT]Z\ 



Note that the eccentricity damping timescale of the sec- 
ond body in the equation above is r e2 . Depending on Q, 
this timescale can appear to greatly exceed r ei . How- 
ever, it is important to keep in mind that in reality, 
variations in fa and fa are coupled because both give 
rise to changes in the planetary semi-major axes. This 
means that tidal dissipation of the inner planet's eccen- 
tricity also damps the outer planet's eccentricity reso- 
nantly. Furthermore, the first and the second planet 
are also coupled through a secular term of the form 
"Hscc oc eie2Cos(n7i — W2), that we have neglected in 
the Hamiltonian. Through this secular interaction, tidal 
damping on e% is translated to e 2 as well (albeit on a 
longer timescale), e ven if there is no d i rect d a mping on 
e 2 (i.e . t £2 = 00; |Wu fc Goldreich| (|2002[); |Mardling 



(20071) 
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Because the dissipation is applied directly on the ac- 
tions, Hamiltonian properties of the solution such as the 
conservation of phase-space area bounded by the orbit 
are destroyed. On a timescale of a ~ few r e , the sec- 
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Since the equations of motion remain linear in z, they 



ond terms in the solutions ( 29 ) will decay away, making 
the phase-space area bounded by the orbit tend to zero. 
This has a number of important physical implications. 
First of all, this removes the dependence of the long term 
(t t z ) solution on the initial conditions. Second, the 
fact that explicit time-dependence of the solution is also 
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lost, suggests that the eccentricity dynamics falls onto a 
fixed point attractor, characteriz ed by constant actions 
(i.e. eccentricities) and angles (Batygin & Morbidelli 
20111. Specifically, assuming that l/r e <C (\a/r)\,\p/r)\) 
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where the involved quantities are given in terms of Ke- 
plerian orbital elements by equations (4), (19) and (20). 
Mathematically, A(j> ps tt arises from the fact that for all 

first-order resonances, /res < 0, while f^s > 0. A physi- 
cal consequence of this fact is that all stationary resonant 
planetary pairs will be apsidally anti-aligned. 

The above solution diverges as rj —¥ and gives positive 
values of e\, e 2 only if 77 < 0. This is because the stable 
equilibrium points of the resonance are always character- 
ized by period ratios ri\ jn% that are larger than the exact 
resonant value. This is a well-know n fact for first or der 
resonances (see for example Ch.9 of |Morbidelli| ( 2002 ) rj 

The solution (29) also illustrates that, beyond the tran- 
sient equilibration period, the eccentricity ratio remains 
constant for all time, since a and /3 are strictly con- 
stant, while the actual eccentricity values depend only 
on 77, i.e. on the proximity of the planets to exact reso- 
nanc^j Because we have restricted ourselves to only a 
linear treatment of eccentricity, this solution fails close to 
exact resonance, where equilibrium eccentricities can be 
quite large. However, this limitation only proves prob- 
lematic in a rather narrow region of parameter space. 

Thus far, we have only considered the relatively fast 
equilibration of orbital eccentricities and critical angles. 
Let us now turn our attention to the truly long-term 
evolution of the system and the associated change in the 
semi-major axes. There are two effects of importance. 
The simpler of the two effects is direct ti dal damping 
of semi-ma jor axes. To leading order in e (Goldreich & 
Soter|p66l ), 

— = -2e 2 -. (31) 

dt J tide T e 

Recall that the eccentricities converge onto quasi-fixed 
points. Thus, in terms of Poincare variables, the tidal 
decay of semi-major axes can be written as: 
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For similar physical planetary parameters (including 
quality factors) and eccentricities, tidal evolution will 
cause orbits to diverge, since r e2 /r ei ~ (k/k — l) 10 / 3 , 

1 This is true only for small to moderate eccentricity values. 

2 Note that at the level of approximation which we have em- 
ployed, the eccentric contribution to \& can be neglected, since 
<E> oc e . Thus, in the definition of r\ in (18) it can be safely as- 
sumed that \P ~ A. 



although both semi major axes drift in the same direc- 
tion (i.e. decay towards the central star). 

The second, more subtle effect is the resonant diver- 
gence of the orbits, forced by eccentricity damping. As 
shown above, tidal decay of eccentricity causes the criti- 
cal angles to collapse onto stable fixed points. However, 
these fixed points are slightly offset from the the actual 
foci. This offset results in a monotonic drift of the semi 
major axes in opposite directions. To understand this, let 
us return to our original formulation of the Hamiltonian. 
An application of Hamilton's equations to Hamiltonian 
Q, evaluated on e and 4> given in (29), yields: 
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where we have made the small angle approximation: 
sm(<fr) ~ (p. Note that the rate of change of the outer 
semi-major axis is positive definite, while that of the in- 
ner semi-major axis is negative definite. In other words 
eccentricity damping always results in the drift of the 
semi major axes in opposite directions, as anticipated 
above. 

The long-term behavior of the resonance can be un- 
derstood by combining equations (32), (18) and (29), to 
yield an equation of motion^] for rj: 

d v 3([h]i(k - 1) + [h] 2 k)(ka 2 T e2 + (k- l)(3 2 T ei ) 



dt r) 2 T ei T e2 
This equation admits the solution 



(34) 



^=(-l) 2/3 {% 3 



<)t 



ik[h] 2 + (k-l)[h]i 



'ei 'ej 

x(fca 2 r e2 + (fc-l)/? 2 r ei )} 1/3 , 



(35) 

where 770 < is an initial condition, corresponding to 
the initial value of 77 for a resonant equilibrium (which 
needs to be negative as shown in (29)) . Note that the 
solution (34) monotonically decreases in time, leading to 
an increase in the absolute value of 77, i.e. an increase 
in the distance between the semi major axes of the plan- 
ets relative to the Keplerian location of the resonance. 
The same 7 7 oc t 1 ^ 3 dependence was observ ed in the sim- 
ulations of Papaloizou & Terquem (2010). Meanwhile, 
the resonant angles, <fi will maintain a near-null libration 
width leading to quasi-constant eccentricity evolution. 

Figure (4) presents an example of such evolution. In 
the case shown, two equal-mass (mi = m 2 — 10 Af) 
planets are started out in exact 2:1 resonance with 
ai = 0.05 AU, ei = e 2 = 0.01, and randomly chosen 
angles. In this calculation, we have set r ei = t &2 and use 
this dissipation timescale as a unit of time (this is vali- 
dated as a result of the adiabatic nature of the evolution) . 

3 Here, the direct tidal and resonant contributions to the evo- 
lution of the semi-major axes have been combined assuming that 
there are n o in direc t ter ms in the disturbing function i.e. the /3's in 
equations ^32\ and ( |33 [ ) are identical. This is true for all first-order 
resonant arguments, except <j> = 2X2 — Xi — In the exceptional 
case, proper account for the indirect terms must be taken (this is 
done in the calculation shown in Fig. 4). 
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Fig. 5. — The fractional extent of divergence away from nominal resonance, A, as a function of the number of elapsed circularization 
timescales, r. The three panels correspond to the 2:1 (A), 3:2 (B) and 4:3 (C) mean motion resonances. The various plotted curves coincide 
with different mass ratios. Particularly, the blue, red and black curves are representative of mi = 3 X 10 — 6 M, mi = 1.5 X 10 — 5 M and 
mi = 7.5 X 10 _5 M respectively. As labeled in the Figure, for each choice of mi, three choices of r»2 = mi/5, ni2 = mi, m,2 = 5mi 
are plotted, with the higher m2 always corresponding to greater A. For all calculations, we set r ei = i~ e2 = t. Note that some systems 
can reach a fractional deviation from exact resonance of up to ~ 20%, suggesting that dissipative divergence of resonant orbits is a viable 
mechanism for production of planet-pairs that reside significantly outside of nominal resonance. 

tion resonances. In the figure, blue curves correspond to 
mi = 3 x 10 _6 M, red curves to m 1 = 1.5 x 10~ 5 M and 
black curves to m\ — 7.5 x 10 -5 M. For each color-coded 
choice of mi, three choices of m,2 = mi/5, m,2 = mi, 
rri2 — 5mi are plotted, with the higher mi always cor- 
responding to greater A. Note that after tjr > 100, 
the more massive examples presented in Figure (5), can 
reside more than ~ 10% away from nominal resonance. 
This points at the viability of creating the near-resonant 
overpopulation observed in the Kepler sample by the the 
mechanism discussed here. 



As above, each panel shows three separate calculations. 
Blue curves represent solutions obtained by numerically 
integrating the non-linear Hamiltonians (TH) and (|6j) in 
pres ence o f tid al dissipation (parameterizeaby equations 
(27) and (31 1), red curves stem from the fully analyti- 
cal framework presented in this section, while the gray 
curves result from an A^-body simulation, where tidal and 
general relativistic inte ractions are accounted for directly 
( Mardling fc Lin|2 002 [) and integrat ed using the Bulirsch- 
Stoer algorithm ( Press et al.||1992 ). As predicted by the 
theoretical arguments above, alter a few (~ 5) circu- 
larization timescales, the system collapses onto a fixed 
state where the critical angles approach their respective 
foci and the variations in eccentricities damp out. Once 
a quasi-stationary configuration is achieved, the orbits 
slowly diverge while the two resonant angles <pi and <p2 
remain in libration which means, strictly speaking, that 
the resonant configuration is maintained (although the 
separatrix ass ociated with the resonance disappe ars at a 



certain r) - see Delisle et al. (20121; Peale (1986)) 



Importantly, when dissipation is applied to a resonant 
pair, the outer orbit drifts outwards, gaining orbital en- 
ergy. This behavior is in contrast with a naive applica- 
tion of standard tidal theory to the individual planets, 
where both planets are taken to drift inwards and facili- 
tates a faster divergence of the orbits. 

As already mentioned above, the long-term evolution 
of the system is adiabatic: the characteristic timescale 
for significant orbital divergence greatly exceeds the res- 
onant interaction timescale. Conveniently, this fact ren- 
ders orbital divergence to be a scale-free process. In 
other-words, the fractional divergence away from exact 
resonance is not explicitly controlled by the actual semi- 
major axes or masses of the planets but rather by the 
mass-ratios (mi/rri2,m/M) and the number of elapsed 
circularization timescales, tjr. Taking advantage of this, 
we have delineated the fractional extent of orbital diver- 
gence, 



Ak:k-1 



n 1 /n 2 — k/(k— 1) 
k/k-1 



(36) 



as a function of elapsed dimensionless time, i/r, for 
an array for planetary mass ratios. These results are 
demonstrated in Figure (5) where the three panels cor- 
respond to the 2:1 (A), 3:2 (B) and 4:3 (C) mean mo- 



4. DISSIPATIVE DYNAMICS OF MULTI-RESONANT 
PLANETARY SYSTEMS 

There is considerable motivation to extend the above 
analysis to systems made of more than 2 planets, where 
each body is in resonance with all of its neighbors, as such 
systems appear to be common in nature. Perhaps the 
best-studied example of a multi-resonant system is the 
Galilean satellites, where both satellite pairs are locked 
in 2:1 mean motion resonances, leading to the libration 
of the Laplace argument. In the collection of confirmed 
extrasolar planets, examples of multi-resonant systems 
include the GL876 sy stem - where the La place resonance 
is directly observe d ( |Rivera et aT1|2010) ), the HD40307 
( Mayor et al. 2009 ) system - which contains three planets 
that reside suspiciously close to a 4:2:1 period commen- 
surability, as well as a few examples in the Kepler data 
set. Furthermore, it has been shown that multi-resonant 
states can serve as good candidates for the initial condi 
tion of the sola r system ( |Morbidelli et al.||2007 Batygin 
fc Brown|[2010| . 

In this section, we shall extend our analytical theory 
of the long-term dissipative evolution of resonant config- 
urations to systems that comprise more than 2 planets. 
As will be shown below, the dynamics of multi-resonant 
systems can be quite rich in diversity, so for simplicity, 
we shall work with a system consisting of three planets, 
keeping in mind that extension to a larger number of 
resonant objects can be accomplished. 

As above, let us begin by writing out the full Hamilto- 
nian. The Keplerian part reads: 



u 



kep 



G 2 M 2 r, 
2Af 



G 2 M 2 m\ 
2A| 



G 2 M 2 m\ 
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while the resonant contribution is: 



G 2 Mm 1 m\ 



res 



i 2T 1 
[Ah 



cos(£i) 



+ fr\ 



(2,in) , 



'2T2 
[A]s 



cos(ei n )) 



G 2 Mm 2 m| (lout) 

w res 



[A]| 



(2, out) 




cos(^ 2 out ) 



2Ta 
[A] s 



cos(^)), 



(38) 



where the four harmonics are: 



! =fc i "A 2 -(fc in -l)A 1 + 7 i 
^ n = fc in A 2 - (fc in - l)Ax + 72 
•«t = A: out A3-(fc out -l)A 2 +7 2 
3 =£: out A 3 -(fc out -l)A 2 +73 



(39) 



and the superscripts "in" and "out" refer to the reso- 
nances of the inner and outer pair of planets respectively. 
Before proceeding further, we note an important differ- 
ence with the formalism developed in the previous sec- 
tion. In the two planet case, dissipation caused both 
critical angles to collapse onto their respective foci. Let 
us examine if similar behavior is possible in the three 
planet case. 

Suppose all four critical angles have evolved to a 
state where d£/dt — 0. In this case, simultaneous 
zero-amplitude libration of d^i/dt — d^ 2 n /dt = and 
d^ 2 ut /dt—d£_ 3 /dt = implies that the apses of the system 
are locked i.e. djx/dt = dj 2 /dt = dj^/dt = d^ sys /dt. At 
the same time, expressing the mean longitude as dX/dt = 
n — dj S y S /dt, the relationship dcj) 2 n /dt — d(f> 2 ut jdt = im- 
plies a strict correspondence among the semi-major axes: 
-k out n 3 + {k in + k out -l)n 2 -(k in -l) ni =0. A configu- 
ration that obeys this relationship is in (or close to) nom- 
inal resonance (e.g. the Galilean satellites). This means 
that away from nominal resonance, only three out of four 
critical angles can reside at their respective foci, while the 
remaining angle will circulate with the frequency 

d£ ciic /dt = -fc out n 3 + (fc in +A: out -l)n 2 -(fc in -l)n 1 . (40) 

Naturally, if the system is far from nominal resonance, 
this circulation is comparatively fast, allowing us to 
drop (i.e. average over) the quic kly varying harmonic 
and reduce the Hamiltonian ( 38 ) to a form that only 
contains three terms. This would further let us con- 
struct new action-angle coordinates, ensuring that the 
momenta conjugated to the three mean longitudes be- 
come constants of motion. However, identifying the cir- 
culating angle is not trivial a-priori, since the calculation 
inevitably depends on the planetary physical parameters, 
and in some cases can have non-linear dependence on ini- 
tial conditions. Thus, unlike the two-planet problem de- 
scribed above, multi-resonant systems should be treated 
on a more case-by-case basis, as the construction of a 
suitable analytical theory for the long-term evolution de- 
pends on the properties of the system. Fortunately, as 
we already showed above, the timescale for the system to 



reach a quasi-stationary state is not much greater than 
the circularization timescale. So the initial transient pe- 
riod of system equilibration can be calculated numeri- 
cally at a mild computational cost. 

Due to the individual attention that multi-resonant 
planetary systems deserve, we shall leave the in-depth 
analysis of detected objects to follow-up papers and in- 
stead limit ourselves to an illustrative example of the 
long-term dynamical evolution of an equal-mass (mi = 
m 2 = TO3 = 10~ 4 Af) planetary system in a 4:2:1 res- 
onance. The aim of the calculation is largely to high- 
light the subtle differences between the evolution of a 
multi-resonant system and the results obtained for a sin- 
gle planetary pair in the previous sections. 

With foresight, we begin with the construction of new 
canonically conjugated coordinates using the following 
generating function (intended for the system at hand): 



F 2 = Ai^i + A 2 * 2 + A 2 * 3 + (fc in A 2 - (k in - l)Aj + 7 i)$i 

(41) 



(r ut A 3 -(A ; out ~l)A 2 +72)<i>2 



+ (k out \ 3 -(k mt ~l)\ 2+l3 )<S> 3 , 
which yields the variables 



*1 


= Ax 


+ {k m - i)$i 






4>i 


= Ai 


* 2 


= A 2 


- fc in $x + (fc° Ut - 


-1)($ 2 4 


-$3) 


"02 


= A 2 




= A 3 


-fc° Ut ($ 2 +$ 3 ) 






^3 


= A 3 




= Tj 


4>i = k in x 2 - 


(fc in - 1 


Ai + 


71 






= Tj 


02 = fc out A 3 - 


- (fc out - 


1)A 2 


+ 72 




$ 3 


= r 2 


03 = fc out A 3 - 


- (fc out - 


1)A 2 


+ 73- 


(42) 



This choice of variables is appropriate when the angle 
£ 2 n is in circulation. Dropping this harmonic from the 
Hamiltonian renders (\&x> ^2, ^3) constants of motion (if 
instead, the circulating angle had been £ 2 ut , the choice of 
\l/ 2 and 4>2 would have been made as in (13), identifying 
k in (13) with k ln , and the angle £ 2 ut would have been 
dropped from the Hamiltonian). 

After some manipulation (as described in the previous 
sections), the Hamiltonian takes on a simple form: 



•H = r; in $i +77 out $ 2 
+ a out V2$2"cos 
where as before, 




Tf" - 
put . 



3([/l]l(fc ln -l)*l-[/l] 2 fc ill *2) 
3([h] 2 (k mt - l)tf 2 - [h] 3 k°^ 3 



(44) 



are the proximities to e xact resonance. The coefficient 
a m is given by equation ( 20 1 and analogously, 



(3° 



G 2 Mm 2 m\ f& out) 

[Mi 



13 V[A] 2 
G 2 Mm 2 ml / r ^ out) 



(45) 



As shown in the previous section, under dissipation the 
system will approach a quasi-stationary state. Once such 
a state is achieved, the corresponding fixed-point orbital 
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Fig. 6. — Dissipative evolution of an equal-mass (mi = m,2 = m-j = 10~ 4 Af) planetary system in a 4:2:1 multi-resonant state over 
t/r = 100 circularization timcscales. As in the two-planet case, the system settles onto a quasi-stationary state. However, the associated 
timescale is somewhat longer: t ~ 10t. As discussed in the main text, only three of four critical angles can equilibrate, while the remaining 
angle is forced to circulate when far from nominal resonance. For the particular setup considered, as shown in panels C and F, the angles 
that tend to their respective foci are £i, £g ut and £3. Meanwhile, the gray (N-body) and green (semi-analytical) dots in panel C show the 
rapid circulation of Panels A, D and E show the evolution of the planetary semi-major axes. In contrast to the two-planet calculation, 
here the drift of d2 is inward rather than outward. Finally, the eccentricity evolution is shown in panel B. Although e\ and eo, settle onto 
quasi-stationary values, is significantly affected by the circulation of £ 2 n j never allowing the eccentricity to fully equilibrate. As before, 
the red curves were obt aine d anal ytic ally, while he blue curves were obtained by numerically integrating the equations of m otio n tha t arise 
from the Hamiltonians | |37[ l and | |38| l, augmented with a simple parameterization of tidal dissipation (i.e. equations | |27|l and | |31| l). The 
gray curves were computed numerically with a direct N-body simulation where dissipation has been taken into account. Note the excellent 
quantitative agreement between the theory and the numerics. 



parameters take on a familiar form: 




1 



resonant drift of the semi-major axes: 
fc in (a in ) 2 



1 

mOUt/7- 

'/ 'e 2 



(46) 



dAi 
~dT 

dA 2 
dt 



dt 



1 



Tex (V in ) 2 

1 (a out ) 2 

n out \2 



(1 - fc out ) 

1 (/3 out ) 2 



r e2 (?7° ut ) 2 t 63 (?7° ut ) 

1 (a out ) 2 1 (/3 out ) 2 
tZ (»7 out ) 2 + tZ iv out ) 2 



(47) 



It is important to recall that we have dropped a quickly 
varying resonant term from the Hamiltonian when de- 
riving these equationfj^] While the dropped harmonic 
will have little long-lasting effect, it will act to introduce 
high-frequency "noise" into the solution, whose ampli- 
tude depends on the proximity of the system to exact 
three-body resonance. Thus, the equilibrium eccentrici- 
ties and critical angles derived here are representative of 
average values. 

Thus far, the behavior inferred from the above equa- 
tions appears quite similar to the case of a single reso- 
nant pair described in the previous sections. However, 
an important difference surfaces when we consider the 

4 Had the quickly varying harmonic been g° ut instead of fo" 1 ) the 



As in the two planet case, the drifts of the innermost and 
outermost planets are inward and outward respectively. 
The migration direction of the second planet, however, 
depends on the relative strengths of the inner and outer 
resonances, since the first term is positive definite while 
the second term is negative definite. Indeed, one could 
envision a set of system parameters (e.g. to 3 -C m^TOi) 
where tidal dissipation leads to a divergence away from 
one set of resonances (increasing |?7 m |) and convergence 
onto another set of resonances (decreasing |?7 out |). In the 
context of such a scenario, conservation of the null phase- 
space area occupied by a quasi-stationary orbit will lead 
to ecc entr icity growth (this can be inferred from equa- 
tions (46)). At the same time, it is important to re- 



coefficients in front of terms containing <I>2 in |43| would have been 
/3 m . Equations U46B would then be modified accordingly. 



call that the presented equations were derived as an ex- 
pansion around nominal resonance location (which is as- 
sumed constant) and thus require dissipation in order 
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to give rise to the corresponding drift of the semi-major 
axes. That is, one could in principle envision a s cenario 
where only r ei is finite, for which equations (47 1 would 



predict a diverging inner pair and a stationary outer- 
most planet, inconsistent with resonant capture (and 
the associated drift of the nominal resonance location, 
d[A]/dt). However, as already pointed out above, the res- 
onant harmonics are non-linearly coupled. Consequently, 
such a situation is atypical in practice, since dissipation 
on a single planet also results in damping of the other 
planet's eccentricities. 

The application of the developed theory is demon- 
strated in figure (6) . For the particular illustrative setup 
considered here, the angles (</>i, 02, ^3) attain a near-focal 
state within t ~ 10r while the dropped harmonic contin- 
ues its circulation as expected. Although all three eccen- 
tricities decay monotonically as before, there is a clear 
qualitative difference in the behavior of e2 compared to 
that of the two-planet case. In particular, e-i never settles 
onto a fixed point, and is instead continuously driven by 
the circulation of which contains 72, an angle conju- 
gated to r2 oc e|. Perhaps unsurprisingly, e\ and are 
not strongly affected by this circulation. 

A more important distinction between the 2-planet 
and 3-planet evolutions is the direction of the second 
planet's drift. Namely, the combined effect of tidal dis- 
sipation and resonant interactions is now to drive the 
middle planet inward, whereas the evolution of 02 was 
positive definite in the 2-planet case. All of this hints at 
the wide variety of possible outcomes and the dynamical 
richness of the multi-resonant interactions in presence of 
dissipative forces. 

5. DISCUSSION 

The primary aim of this work has been to formulate a 
simple, physically intuitive analytical theory for the dis- 
sipative divergence of resonant orbits. We began with a 
purely conservative treatment of a single resonant pair 
and showed that at sufficiently low eccentricities and 
limited libration amplitudes, resonant dynamics can be 
treated with a linear, integrable approximation to the full 
resonant Hamiltonian. We then introduced simply pa- 
rameterized tidal dissipation into the equations of motion 
and showed that the system tends to a quasi-stationary 
state over a few eccentricity circularization timescales. 
The collapse of the critical angles onto near-focal values 
in turn results in a divergent drift of the semi-major axes 
such that the outer orbit continually gains orbital energy 
while the inner planet's orbit decays. We subsequently 
showed how the developed formalism can be extended to 
multi-resonant systems. However, we have limited our- 
selves to a single illustrative example of the evolution 
of a system near a Laplace-like resonance, as we argued 



that the parameter space available to multi-resonant sys- 
tems is quite large, rendering individual modeling more 
cost-effective. 

Overall, our results point at the distinct possibility 
that the dynamical architectures of numerous detected 
systems, whose orbits seem to lie outside of resonance on 
the basis of the observed orbital periods, are a resu lt of 
resonantly-aid ed dissipative divergence of the orbits ( Pa- 
paloizou 2011 ), and thus comprise a number of important 
implications. First, the explanation we propose suggests 
that protoplanetary disks are indeed conducive to form- 
ing resonan t planetary systems, whose long-term survival 
is assured (Cresswell & Nelson 2008[ ). In combination 
with precise quantitative modeling, this constraint can 
likely yield important new insights into understanding 
the physical structure and evolution of protoplanetary 
disks (e.g. weakly turbulent). 

Second, as shown in section 3, depending on the mass 
ratio and the elapsed time, resonant orbits can evolve 
up to tens of percent away from nominal resonance. If 
such extreme evolution is common, it is possible that 
many planetary systems are actually in resonance even 
if their orbital periods are apparently not in commensu- 
rability. In particular, we expect the period-ratio statis- 
tics of newly-formed planetary systems to cluster more 
clearly arou nd resonant values tha n those of an evolved 
sample (see Fabrycky et al. (2012) for an in-depth dis- 
cussion of the current data setf. 

Third, the fact that the time-dependence of the 
orbital divergence is related to the tidal circularization 
timescale can be used to infer from the observed period 
ratio how many circularization timescales a given 
system has evolved through, if the age of the system is 
known. Such information is vital for constraining un- 
observable parameters of extra-solar plane tary systems 
such as the planetary tidal quality factor (Goldreich & 
Soter 1966), whose origin remains largely unexplained 
and is among the most poorly constrained values in 
astrophysics. Although the above arguments hinge 
on the observationally elusive characterization of the 
physical planetary properties, we can certainly expect 
the data to improve continuously over the coming years 
allowing for these calculations to be executed, eventually. 
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